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Abstract. We describe an explicit in time, finite-difference code designed to simulate black 
holes by using the excision method. The code is based upon the harmonic formulation of 
the Einstein equations and incorporates several features regarding the well-posedness and 
numerical stability of the initial-boundary problem for the quasilinear wave equation. After 
a discussion of the equations solved and of the techniques employed, we present a series of 
testbeds carried out to validate the code. Such tests range from the evolution of isolated black 
holes to the head-on collision of two black holes and then to a binary black hole inspiral and 
merger. Besides assessing the accuracy of the code, the inspiral and merger test has revealed 
that the marginally trapped surfaces contained within the common apparent horizon of the 
merged black hole can touch and even intersect. This novel feature in the dynamics of the 
marginally trapped surfaces is unexpected but consistent with theorems on the properties of 
these surfaces. 



PACS numbers: 04.25.Dm, 02.70.-c, 02.70.Bf, 02.60.Lj 

1. Introduction 

The numerical calculation of the inspiral and merger of a binary black hole system attracted 
early attention because its dynamic and strongly relativistic gravitational fields were expected 
to play a major role in astrophysics and to provide an excellent arena for studying the 
fully nonlinear behavior of gravitation. It has long been known £0 that the final decay, 
coalescence, and ringdown of such a system is a very strong source of gravitational radiation. 
Because of their complexity and nonlinearity, and the lack of any continuous symmetries, the 
Einstein equations cannot be solved analytically for these systems. Rather, direct numerical 
integrations must be used. 

Building on an earlier attempt by Hahn and Lindquist (2j, the first successful numerical 
simulations of binary black hole systems were performed in the 1970s by Cadez 
Eppley |4), DeWitt, and their colleagues |6] |7). These simulations were restricted 
to axisymmetry, and used the Arnowitt-Deser-Misner (ADM) formulation of the Einstein 
equations |8|. However, the ADM equations are now known to be only weakly hyperbolic [9|, 
so they are not suitable for long-term numerical evolutions despite large-scale efforts ifTOl . 

Starting with the first long-term evolutions by Pretorius ifTTI . there has been remarkable 
progress in the simulation of binary black holes. Several groups have since tracked the 
inspiral and merger phases using codes solving a conformal-traceless formulation of the 
Einstein equations and treating the black holes either by the excision method [121 or by 
the moving-punctures method lfl3l [141 [T5l \W\ . Since Pretorius' original work, there has 
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also been substantial progress in developing mathematical theorems which establish the 
well-posedness of the harmonic initial-boundary value problem and the stability of its finite 
difference approximations. We have incorporated some of this theory in developing an explicit 
in time, finite-difference harmonic code to treat the excision problem. We present here some 
promising initial results indicating that the binary black hole merger can be treated using 
excision with a minimal amount of dissipation and with sufficient accuracy to reveal new and 
interesting dynamics of the individual marginally outer trapped surfaces (MOTS) contained 
within the common apparent horizon. 

Harmonic coordinates were first introduced by deDonder ifTTl to reduce Einstein's 
equations to 10 quasilinear wave equations and they were later extensively developed by 
Fock [18] and used by Choquet-Bruhat 1 19] to give the first well-posed version of the Cauchy 
problem for the gravitational field. In a harmonic gauge the spacetime coordinates (viewed 
as a set of 4 scalar functions) satisfy the curved-space wave equation V ' c V c x^ = 0. The 
principle part of the Einstein equations then reduces to a second-order hyperbolic form or 
to a first-order symmetric hyperbolic form, for which there is an extensive mathematical 
and computational literature. Many researchers have since implemented numerical evolution 
schemes for harmonic formulations of Einstein's equations J20]|2T][22 23 24l l25ll26lfTTll27l 
[28] |29] [30] and the related Z 4 formulation l3Tll32l . 

The AEI harmonic code presented here has its roots in the Abigel code l22l |23~1 l25l . 
a second-order accurate, finite-difference code which incorporates theorems establishing the 
well-posedness and numerical stability of the harmonic initial-boundary value problem. In 
addition, the AEI code incorporates a black-hole excision algorithm which allows for motion 
of the excised region across the grid, a superluminal evolution algorithm and, except for 
regions near the boundaries and for some aspects of the mesh refinement, a fourth-order 
accurate finite-difference approximation. It also utilizes an apparent horizon finder l33ll34l . 
vertex-centered mesh-refinement techniques 11351 . black hole initial data sets and other 
features of the Cactus computational toolkit 11361 necessary for black hole simulations. A 
Runge-Kutta integration is used to carry out an explicit time evolution. Note that this differs 
from the approach of Pretorius who uses a pointwise Newton-Gauss-Seidel relaxation scheme. 

Although some aspects of the code are still under development, most notably the excision 
boundary and the outer boundary treatments, we report here a series of testbeds carried out 
to validate the code. Such tests range from the evolution of isolated black holes to the head- 
on collision of two black holes and then to a binary black hole inspiral and merger. While 
these tests have now become standard, we have found a new feature in our study of the 
binary inspiral and merger. In particular, after a common apparent horizon has formed, our 
simulations show that the two individual MOTS approach and continue on to intersect each 
other. This novel feature in the dynamics of the MOTS is unexpected but not unreasonable if 
these surfaces have to maintain their smoothness as they meet. It is also consistent with recent 
theorems concerning the properties of MOTS (L. Andersson, private communication). 

The paper is organized as follows. In Sect. |2] we present the system of evolution 
equations and describe its reduction to first-order in time form, as well as our use of gauge 
conditions and boundary conditions. In Sect. [3] we give a brief discussion of the numerical 
techniques we have employed. Section|4]collects the results of our tests and the calibration of 
the code's accuracy. Finally, Sect. [5]summarizes our results and the prospects for future work. 
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2. Harmonic Evolution System 

2.1. The evolution system 

Although completely redesigned, the AEI harmonic evolution code is based on the work 
presented in Babiuc, Szilagyi et al. l22l l23l l25l with modifications to allow for a smooth 
transition from subluminal to superluminal evolution, together with higher-order finite- 
difference operators and the possibility of excising an arbitrary portion of the grid (moving 
black-hole excision). 

In a generalized harmonic gauge 11371 . the coordinates x p = (t, x l ) = (t, x,y, z, ) satisfy 
- V a V/ = = F* 1 , (1) 

where 

T» := = —±=:d v r u , (2) 



g 

with gauge source functions F tM (x p , g pa ) (which may depend on the spacetime coordinates 
and the metric) and with the densitized 4-metric 

r v ■■= V~gg" v (3) 

playing the role of the basic evolution variable. In this harmonic formulation, the constraints 
reduce to the gauge condition 

C := T p - F p = , (4) 

and the evolution system is based upon the reduced Einstein tensor 

Ei> ,v ._ Qiiv _ V (p r ^) + _ g ^ VaT a _ (5) 

Here Y v is treated formally as a vector in constructing the "co variant" derivative VT". When 
the constraints (|4]) are satisfied, this gives rise to a hyperbolic evolution system 

E ?u = _ v (m f ^) + IgWVpFP. (6) 

Provided the gauge source functions do not depend upon derivatives of the metric, they do 
not enter the principle part of the system and do not affect its well-posedness or numerical 
stability. The evolution system © takes the specific form of quasilinear wave equations 

d P {g pa d a ~ 9 n - 2^g pa g TX ^ pT K x - V=g(d p9 n(d.gn + ^={d p9 n{d«g) 

+ 1<T U^=.{d p g){d a g) + yf^drfT + -^= g {d a g)d p gP^j 

+ 2^g~V^F» ) - y/^gg^VpF" + V^gA^ = , (7) 
where we have included the possibility of a constraint-adjustment term 

#":=C^(/, V ,3rV). (8) 
i.e., a term which vanishes when the constraints are satisfied and which does not affect the 
principle part of the evolution system. 

Note that we do not explicitly enforce the harmonic constraints (j4| during the evolution. 
Instead, we invoke the Bianchi identities which imply wave equations of the homogeneous 
form 

g pa d p d a C p + UydpC + M%C a = , (9) 
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where the matrices L and M are functions of the metric and its first and second derivatives. 

Given constraint-preserving initial and boundary conditions, the uniqueness of the 
solutions to (O guarantees that the harmonic constraints be conserved during the evolution. 
On the other hand, constraint-preserving initial data also requires that the initial Cauchy 
data satisfy the standard Hamiltonian and momentum constraints. Also, since the harmonic 
constraints imply evolution equations for the lapse and shift, the only remaining free initial 
data in addition to the usual Cauchy data (the 3-metric and extrinsic curvature of the Cauchy 
hypersurface) are the initial choices of lapse and shift and of the gauge source functions. Note 
that an initial choice of the gauge source functions is effectively equivalent to a choice in the 
initial evolution of the lapse and shift. 



2.2. Constraint Adjustment and Damping 

The constraint-adjustments implemented in the code are those investigated by Babiuc et al. ll25l 
and have the general form 

AIU , _«L =C p a *»> + a2 ° P ^ p t C»C V - -^=C^V v) t , (10) 

where the a, > are adjustable parameters (in the runs reported here we have set aj = 1), 
e CTT is the natural metric of signature (++++) associated with the Cauchy slicing and e 
is a small positive number chosen to ensure numerical regularity. The effects of these 
adjustments in suppressing long wavelength instabilities in standardized tests for periodic 
boundary conditions have been discussed by Babiuc et al. ||251 . 

In particular, the first and second terms in the adjustments ([Toll have been shown to be 
effective in suppressing constraint- violating nonlinear instabilities in shifted gauge-wave tests. 
The third term in (fT0|) . on the other hand, was first considered in ll38l and leads to constraint 
damping in the linear regime. Although it has been used effectively by Pretorius lfTTll27l in 
black-hole simulations, it was not effective in the nonlinear regime of the shifted gauge-wave 
test |25l . 

We also note that the work reported in ll25l has shown that adjustments which scale 
quadratically with C p (or with higher powers) take effect too late to counter the growth of a 
constraint-violating instability; and this has lead to the specific form for the denominator of 
the second term in ( TTOb . 



2.3. Gauge Conditions 

As noted earlier, the gauge source functions F 1 - 1 may be chosen to be arbitrary functions of the 
spacetime coordinates and metric. They can be viewed as differential gauge conditions on the 
densitized metric. This serves two important purposes. Firstly, it allows for convergence tests 
based upon a known spacetime, whose analytic metric (x p ) is specified in a non-harmonic 
gauge, by choosing 

F» = — T t=d v g% ) . (11) 

Using these analytic gauge source functions, in combination with initial and boundary data 
consistent with the analytic solution, gives rise to numerically evolved spacetimes that are 
identical to the analytic solution up to discretization error. This is how the convergence 
tests reported here have been carried out for the Schwarzschild spacetime expressed in 
(non-harmonic) Kerr-Schild coordinates. Secondly, and most importantly, the gauge source 
functions can be used to avoid gauge pathologies. 
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A major restriction for the form of the gauge source functions is that they cannot depend 
on the derivatives of the metric. In particular, they cannot depend on the location or shape 
of the MOTS and this is a problem when moving black holes are present, and where it 
is important for the coordinates to be able to "respond" to the black hole motion. In our 
simulation of binary black holes, we have used the gauge source function 

F* 1 = ~J= - V*n , (12) 

where rf v is the Minkowski metric and where w = ^(x 1 ) is a smooth, spherically symmetric, 
time-independent weighting function with ui = 1 over most of the computational domain, 
but with to = in some neighborhood of the outer boundary. When spatial derivatives are 
neglected and w = 1, the resulting gauge condition takes the simpler form 

d t g tfi = -(9 tfi - V*") , (13) 

showing that it forces the densitized lapse and shift to relax to their Minkowski values. 

In our first attempts at binary black hole simulations, we have found that this choice 
of gauge source function keeps the lapse and shift under reasonable control. Similar 
choices of gauge source functions have been used with success in other binary black hole 
simulations lfTTll27ll . 



2.4. Boundary Conditions 

Our evolution domain has a timelike outer boundary and a smooth, spacelike excision 
boundary inside each MOTS. The harmonic evolution system, in the second-order form 
(0, consists of quasilinear wave equations whose characteristics are identical to the null 
directions determined by the metric. As a result, all characteristics leave the spacelike excision 
boundaries and no boundary conditions are necessary (or allowed). 

At the timelike outer boundary, any dissipative boundary condition for the wave equation 
with shift leads to a well-posed initial-boundary value problem (IBVP). Such dissipative 
boundary conditions were first worked out in the one-dimensional (ID) case l39l l40l HTl 
and general results for the 3D case have been discussed recently in [23 42 1. For a boundary 
with normal in the +x direction, such dissipative boundary conditions have the form 

[(1 - K)dt + Kg zp d p ] r v = <f\ < k < 1, (14) 

for each component g^, where q^ v is the boundary data. The choice k = gives a Dirichlet 
condition and k = 1 gives a Neumann condition. Dirichlet and Neumann conditions are 
marginally dissipative in the sense that they are purely reflective for modes with q^ v = 0. A 
strictly dissipative Sommerfeld-type condition arises when k is chosen so that the derivative 
in the left hand side of ( TPfl i lies in the outgoing null direction. 

In order for the IBVP to be constraint preserving, the boundary data q^ v must be assigned 
to enforce a homogeneous, dissipative boundary condition on the constraints C 7 '. Then, with 
proper initialization, the uniqueness of solutions to eqs. © ensures that the constraints are 
satisfied throughout the evolution. The first proposal for such constraint-preserving boundary 
conditions for the harmonic system consisted of a combination of 3 Dirichlet and 7 Neumann 
conditions on the components ofg^ ll22l . However, numerical studies l23l showed that these 
Dirichlet-Neumann boundary conditions were effective in carrying the signal off the grid but 
that their marginally dissipative nature reflected the noise and gave poor results in highly 
nonlinear tests. 

The first example of strictly dissipative constraint-preserving boundary conditions which 
would in principle allow numerical error to leave the grid, was given for a tetrad formulation 



An explicit harmonic code for black-hole evolution using excision 



6 



of the Einstein equations by Friedrich and Nagy [43]. Constraint-preserving boundary 
conditions of the Sommerfeld type which lead to a well-posed IBVP for the non-linear 
harmonic formulation have subsequently been formulated l44l . These Sommerfeld-type 
boundary conditions have been incorporated in a numerical code which gives vastly superior 
performance in nonlinear test problems than the Dirichlet-Neumann scheme |24| . However, 
we have not yet implemented these condition here; instead we have used a naive version 
of Sommerfeld boundary conditions which not only is not constraint-preserving but whose 
numerical implementation is only second-order accurate. In addition, as explained further in 
Sect. [3j the code uses summation-by-parts (SBP) difference operators which are fourth-order 
accurate in the interior but only second-order accurate in the vicinity of the outer boundary. 

2.5. Reduction to First-Order in Time 

In contrast to what is done in lfTTll27ll . where the harmonic evolution system is second-order 
in time, we find it convenient to discretize and use the method of lines to time-integrate an 
evolution system which is first-order in time. We note that the reduction to first-order in time 
can be done in a number of ways, some of which may have very different stability properties 
when discretized (see |Appendix A| for a discussion). 
Here, we introduce the auxiliary variables 

CT ■= 9 U dtr v + wg^diF" , (15) 

where w — w(x l ) is a smooth weighting function with w = 1 over most of the computational 
domain but with w = in a neighborhood of the outer boundary. We rewrite the 
definition (|15|) to obtain the time derivatives of g^ v in terms of and spatial derivatives of 

d t r u = ^ (q^ - w %r v ) ■ (i6) 

Next, we use the identity (fTo*|) to re-express the principle part of the harmonic evolution 
equations (JT]) as 

d P (g pa d a rn = dt& v + (i - w)d t (g t %rn + + , hd 

and after using (fT6|) to convert all time derivatives of g^ v in (fT7|) and (jT]) into spatial 
derivatives, we obtain an equation of the form 

dtQT = F^(g, 8^, d^g, Q, t Q) . (18) 

Equations (fTo*)) and (fl8j) represent our basic evolution equations for the field variables g p,v 
and Q^ v , respectively. 

3. Numerical Implementation 

3.1. Finite- difference algorithms 

The code solves the finite-difference equations on a Cartesian grid with finest resolution 
Ax 1 = h, using a cubic outer boundary and with excision boundaries for each black 
hole. Vertex-centered mesh refinement is applied using the Carpet driver J35l . within 
the framework of the Cactus computational toolkit 1361 . The time evolution is carried out 
by the method of lines using a fourth-order Runge-Kutta scheme, with a fifth-order spatial 
prolongation and a second-order time interpolation to provide fine-grid boundary data at 
mesh-refinement boundaries. 
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While the bulk of the code uses fourth-order accurate centered difference operators to 
approximate spatial derivatives, in a neighborhood of the outer boundary we approximate 
the spatial derivatives by diagonal norm SBP difference operators of fourth-order interior 
accuracy and of second-order accuracy at the boundary. More specifically, on a grid xj = 
xq + ih with boundary at xq, these operators, as described by Mattsson and Nordstrom 1451 . 
are 



~ Q/m - , d9) 

1 / 4 59 59 4 \ 

h ri3 /[4] + 86 /[3] _ 86 /[1] + 43 /[0] ) ' (20) 



(d x f)i=i 

(0*/)<=8 \ (-^/[5] + §/[4] - jjg/[2] + |j/[0]) , (21) 



and 



(^/) i= i^^(/ [2] -2/ [1]+ / [0] ) , (22) 

, s2n 1 ( 4 59 110 59 4 \ 

[p x t)i=2 -> ^ I — + 43 /[3] - -43- J[2] + 43 /[i] - ^/[o] I > (23) 

,2 - 1 / 4 „ 64 „ 118 „ 59 



(0 x /)i=3 ^ + 49^] - 49/p] + 49/p] " 49/M J ■ (24) 

For grid points at the outer boundary, on the other hand, all components of Q^ v are 
updated using a flat-spacetime, homogeneous Sommerfeld boundary condition. With the 
outer boundary located in the weak-field regime, a Sommerfeld condition applied to Q^ v 
is equivalent to setting the Sommerfeld derivative of g^ v to the value of this Sommerfeld 
derivative at t = 0, as determined by the initial data. (By implication, our boundary algorithm 
is compatible with the initial data.) This boundary condition is effective in maintaining 
numerical stability. However, it is not constraint-preserving and is a prime target for future 
code improvement. 

In the part of the computational domain near the outer boundary where w = 0, the 
evolution algorithm for Q and g reduces to a fourth-order version of the subluminal evolution 
algorithm for evolving the wave equation with shift discussed in |42j . This algorithm is known 
to be unstable in the region where the shift is superluminal (e.g. near the excision boundary). 
In this superluminal region, we set w = 1 so that the definition of Q [cf. eq. Q3) makes Q 
the derivative of g in the normal direction to the Cauchy hypersurfaces, which stabilizes the 
algorithm. A more detailed discussion of this subluminal-superluminal blending is discussed 
in I Appendix A| 

Another important ingredient of our code is numerical dissipation. We find this essential 
in keeping the algorithm stable in the neighborhood of the excision domain. In addition, this is 
also helpful in killing off high-frequency noise generated at the mesh-refinement boundaries. 
In the interior of the grid, numerical dissipation is added at 0(h 5 ) in the form 

d 2 t r v - + ^h 5 ]T e t {D +l D^fd t r v , (25) 

i 

where D±i are the forward and backward difference operators in the x % direction and ej is 
a smooth weighting function. In the neighborhood of a face of the outer boundary with 
normal in the .T-direction, we set e x = so that the dissipation applies only in the tangential 
directions. In carrying out convergence tests for a Schwarzschild black hole, we choose 
e, = 0.2 outside the apparent horizon (AH) and e; = 2 inside the AH (except for a transition 



An explicit harmonic code for black-hole evolution using excision 



8 



region). In the two black hole simulations, we choose ej = 1 outside the AH and ei = 2 inside 
the AH. 

3.2. Moving Excision 

The excision algorithm is driven by the apparent horizon finder algorithm ll33l[34l . Strictly 
speaking, this algorithm searches for MOTS, regardless of whether these are apparent 
horizons or not. We recall that MOTS are defined as smooth, compact 2-dimensional surfaces 
whose outgoing normal null geodesies have zero expansion. With respect to a 3 + 1 foliation, 
the apparent horizon is defined as the 3-dimensional hypersurface traced out by the outer 
boundary of the trapped region in each time slice. If sufficiently smooth, the apparent horizon 
is foliated by MOTS. 

A smooth spacelike boundary is used to excise a region inside each MOTS, resulting 
in a jagged boundary in the Cartesian grid. The excision boundary is chosen to be centered 
inside the MOTS and scaled in coordinate size to be 0.7 the size of the MOTS in the binary 
simulations and 0.8 in the Schwarzschild black hole tests. 

We keep the same interior evolution stencil near the excision boundary by introducing 
the necessary ghost points. Because the dissipation operator d25l l would require an excessive 
number of ghost points we replace it with a third-order form /i 3 (D + ;Z?_;) 2 near the excision 
boundary. Values at the required ghost points are supplied by an extrapolation scheme which 
was proved to be stable for the case of a boundary aligned with the grid 0421 . We have 
generalized this to the case of a generic smooth boundary in a Cartesian grid following 
the "embedded-boundary" method developed by Kreiss and Petersson ll46l for formulating 
a stable Neumann condition. More specifically, we construct a vector v % by taking the flat- 
space displacement from the centroid of the excised region to the current position, and require 
that 

X>*Abi)V = and Yj^D^fCt v = , (26) 

i i 

where the one-sided differences D±i correspond to the sign of v % . The extrapolation condition 
d26l i is applied iteratively at the points near the boundary until the full stencil of ghost points 
is updated. 

Extra care needs to be paid to the identification of ghost points when the excision 
domain is moving across the grid. In particular, grid-points that become interior points at 
t := to + NAt but were ghost points at t^ -1 need to be treated as ghost points, i.e., the 
excision algorithm must fill them with extrapolated values during the time-integration from 
t N ~ x to t . However, when evolving from t N to t N+1 , these same grid-points need no longer 
be treated as ghost points and can then be labeled as evolution points. 

4. Code Tests 

We now present a series of tests to validate our code's stability and accuracy, starting with 
the evolution of Schwarzschild black holes using static gauge source functions, followed by 
the head-on collision of two equal-mass black holes and finally the inspiral and merger of the 
QC-0 initial data set for binary black holes. 

4.1. Single Black Holes 

We have evolved Schwarzschild black holes for long periods of time to test the stability and 
accuracy of our evolution algorithm, excision scheme and outer boundary treatment. To study 
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Figure 1. Error r AH /M — 2 in the areal radius of the apparent horizon as a function of time 
for evolutions of Schwarzschild initial data at 4 different resolutions. The main figure shows 
r AH /M — 2, while the inset shows (O.IM//1) 4 X (r AH /M — 2), demonstrating that the error 
shows fourth-order convergence to zero as the resolution is increased. 
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Figure 2. Scaled harmonic constraint (O.lM/ft) 4 X C° at t = 200M along the x axis, for 
evolutions of Schwarzschild initial data at 4 different resolutions. 
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Figure 3. Masked 2-norm of the harmonic constraint | |C° 1 1 a as a function of the finest grid 
resolution h, at selected times, for evolutions of Schwarzschild initial data at 4 different 
resolutions. Note that the t = 200M, t = 500M, and t = 1000A/ curves are almost 
coincident (because the norms are almost time-independent over this range of times). 

the code's convergence with resolution, we made evolutions with the finest grid spacing set to 
h = 0.100M, 0.080M, 0.0666M, and 0.050M. Each evolution used a similar grid structure, 
with 5 levels of nested 2:1 mesh refinement with the finest grid (spacing h) extending from 
the origin to 4M + 9 h, the next coarsest (spacing 2 h) from the origin to 8M + 9 h and so 
on, up to the coarsest grid (spacing 16 h) from the origin to 64M. The tests were carried 
out in octant symmetry. For representative runs, we cross-checked that identical output was 
produced in full-space simulations. 

Each evolution was run for 1000M, with no sign of instability for this full duration. We 
measure the accuracy of these evolutions by monitoring the apparent horizon's areal radius 
and the harmonic constraint C°. (Convergence results for the error in <? 00 are very similar to 
those for C°.) 

Figure [TJ shows the error in the apparent horizon areal radius r AH /M — 2 for these 
evolutions. After an initial transient, each evolution displays a slow growth (roughly linear in 
time) in the apparent horizon. The growth rates are small, and show fourth-order convergence 
to zero with the resolution. Figure [2] on the other hand, shows the magnitude of C° along 
the x axis at t = 200A/. It is clear that while there are regions of excellent pointwise fourth- 
order convergence, there are also regions where mesh-refinement and other effects make the 
convergence more problematic. 

To gain a clearer picture of the overall convergence behavior of these evolutions and 
because the error is by far and large dominated by a small set of grid points near the excision 
boundary, we use a "masked" 2-norm | j C° | j 2 defined as the root-mean-square norm of C° 
over only those grid points which are outside the apparent horizon and do not have any finer- 
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refinement-level grid points overlaying them. For example, all of the finest grid outside the 
apparent horizon would be included in the masked norm, but only that part of the grid with 
spacing 2 h which is outside 4M +9 h would be included, and similarly for the other coarser 
refinement levels. 

After a short initial transient (lasting a time comparable to the outer-boundary crossing 
time), 1 1 C° 1 1 2 is essentially time-independent in each run. Figure [3] shows the variation of 
1 1 C° | ! 2 with resolution at selected times. As can be seen, the norm decreases with resolution, 
with a convergence exponent which is always > 3, but varies with resolution. We attribute 
this behavior to not having sufficient resolution to see the asymptotic convergence behavior 
as introduced by the dissipation operator. 

Similar results have also been obtained for spinning black holes, which we have evolved 
with very high spins. Simulations of Kerr black-holes with spins J/M 2 up to 0.99 have shown 
no apparent signs of instability up to t = 100 M. We have not yet investigated them further. 

4.2. Binary Black Hole Mergers 

4.2.1. Head-on collision The binary black hole collision we have carried out is an equal- 
mass, non-spinning head-on collision of two Brill-Lindquist black holes of masses 0.5 each, 
located on the z axis at z = ±1.16. 

The test was carried out in octant symmetry, using a pure harmonic gauge {i.e., = 0). 
The numerical grid was set up using nine levels of refinement, with grid-step hS n > = 
3.2 x 2~ n , n = 0, . . . , 8 giving a grid-step of h = 0.0125 on the finest refinement level 
and an outer boundary at 14AM. The chosen grid setup is such that the initial black holes are 
contained within the bounding box of the finest grid and hence no motion of the finer grids is 
needed. 

After the merger, when the final horizon has a an ellipsoidal shape with maximum and 
minimum radii in a ratio r m i n /r max > 0.6, the finest refinement level is dropped. This is 
justified by the fact that the coordinate-radius of the final apparent horizon radius is more than 
double that of the two individual MOTS. Overall, the simulation shows no sign of instability 
and Fig. |4] shows the £ = 2, m = 2, even-parity multipole of the Zerilli function Q^ as 
extracted at R = 60M, indicating a well captured quasi-normal ringing of the black hole. 

4.2.2. QC-0 inspiral We next consider the evolution of the QC-0 initial data defined by l47l 
and constructed with the highly-accurate spectral elliptic solver of 1481 . In contrast with what 
was done for the head-on collision, we here use the gauge source function (fT2l in order to 
keep the lapse from collapsing and the shift from forming large coordinate distortions in the 
strong-field region. While helping in terms of stability, these gauge source functions also lead 
to a coordinate-growth of the MOTS. 

The simulation started out with nine levels of mesh refinements, with the spacing on 
the individual refinement levels being h^ n ' = 2.048 x 2~", n = 0, . . . , 8 and yielding a 
resolution of 0.008-A/ on the finest refinement level, with an outer boundary at « 160A/. The 
simulation was carried out in the x > 0, z > quadrant, taking advantage of the reflection 
symmetry across z = and the ir/2 rotation symmetry around the z axis. Around each black 
hole we have used a set of refinement levels of size i'"' = R AH (t = 0) X 2^ 9- ™^ so that the 
finest grid (n = 8) would initially be twice the size of the apparent horizon. This refinement 
structure follows the motion of the black holes across the grid and, as the coordinate growth 
of the MOTS takes place, we adjust our grid-structure by first dropping the finest and later the 
second finest refinement level. After dropping the second finest level, we re-adjust our grid- 
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Figure 4. Zerilli waveform for the head-on problem, extracted at R = 60M, The outer 
boundary in this test was at L = 14AM. 

structure by setting = i? AH (<) x 2 { - 7 - n \n = 0, . . . , 6, keeping the ratio L^/R An (t) 
constant as the MOTS grows further. 

At the time of the formation of the common apparent horizon the finest grid has 

= 0.032M. The individual MOTS were tracked as long as possible. Eventually the 
horizon finder algorithm hits the excision domain of the other black-hole, at which time the 
individual MOTS are lost. 

Figures [5] and [6] reveal an interesting feature of the MOTS dynamics following the 
formation of a common apparent horizon. The individual MOTS continue to inspiral until 
they touch (apparently at a single point) and then overlap. Note that, by definition, at the 
time of touch (and of later intersection) of the individual MOTS, there already exists a third, 
common MOTS, which is now the apparent horizon. We measure a coordinate velocity 
of v as 0.16 of the MOTS as they approach each other. This phenomena, including the 
approximate value of v, was reproduced using a variety of numerical parameters, including 
varying resolutions, regridding rules, and outer boundary locations. Finer grids allow longer 
tracking by the apparent horizon finder, which in turn leads to larger overlaps between the 
individual horizons. The figures we present here were obtained with our best-resolved run. 

We have measured the convergence of our QC-0 test by comparing results obtained at 
three different resolutions. At the time of the overlap of the MOTS, the grid-steps were 
h h = 0.032A/, h m = 0.0AM and h = 0.044M. The MOTS touching time for the three 
resolutions was as 10.73M, T m as 10.3671/ and T\ = 10.29. As a first measure of the 
numerical error, we checked the convergence of ||C ||oo at t as 10.956. With the excision 
being a major source of numerical error, the norm is a reflection of the constraint error 
near the excised points, including those within the region of overlapping MOTS. Overall, we 
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Figure 5. are (x, y, t). The overlap of the MOTS is clearly visible. Note that these surfaces 
are actually smooth everywhere; their apparent non-smoothness in some parts of this figure is 
an artifact of the perspective projection. 



measure that 



IjgfUc 

\K\l 



2.25 



(27) 



Such a convergence order can be explained by bearing in mind that the derivatives of the 
metric contained in the constraint C° are computed using five-point centered stencils. Near 
the excision boundary, however, two of those five points are updated by extrapolation, using 
the available three points. This effectively implies that the quantity C° is computed near the 
excision boundary using a three-point non-centered, first-order derivative stencil, which has 
an error of 0{h 2 ). 

As an additional measure of the convergence order we have computed the L2 norm of the 
differences between the coordinate-shapes S of the MOTS for the different gridsizes. Given 
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Figure 6. Intersection of the MOTS with the 2 = plane at the time when the individual 
MOTS first touch (t ss 10.73M ) and just before the individual MOTS are lost (t R3 11.47M); 
for the latter time the individual horizon-finder angular grid points are also shown. At the times 
shown here the finest grid has = 0.032Af ; this is illustrated in the plot legend. Notice 
that the MOTSs are well-resolved on both the Cartesian and angular grids at all times, and that 
they clearly overlap at the latter time. 



the fact that the three runs had shape information available at different discrete time values, 
first we interpolated in time, from all three sets of data, to a common set of time-slices 

n — 1 

t n = imin + _ ^ X (i max - t min ), 1% = 1 • • • N. (28) 

The value of N was not critical in this convergence test, given the smooth behavior of the 
surfaces. We used N = 1000 for the individual MOTS and N = 100 for the common 
apparent horizon. 

At each time-level t n , for each resolution, we had a set of grid-points labeled by their 
angular grid-index [/, J] as well as Cartesian coordinates. (Note that for convenience we used 
the same angular grid-resolution of the horizon-finder algorithm for all our runs.) For a given 
choice of t m [ n , i max and N we define the L2 norm of the error in the shape of the MOTS to 
be 



\Si-S h \\ 2 



£[(x|(tn,J, J)-x\(t n ,I, J)f /(NjNjN) .(29) 

\ I,J,n i=l 



We set i max = 11.0 based on the maximum physical time for which horizon data was 
available for the coarsest run. The convergence rate for the individual MOTS (with t m j n = 0) 
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was 



\\S m — Sh\\2 



( 




) 



2.42 



(30) 



The accuracy of the touching and overlapping MOTS was measured by restricting the Li 
norm to t > t m [ n = 10.4. The measured convergence rate is C(ft, 2 59 ), while the common 
(apparent) horizon shape was convergent to 0(h 2 - 57 ). 

Overall the accuracy of our QC-0 test is between second and third-order, which is lower 
than what was measured for the case of an isolated Schwarzschild spacetime. As mentioned 
above, the most likely causes of the larger error are the close proximity of the moving 
excision algorithm to the apparent horizon finder interpolation stencils and the 0(di i ) time- 
interpolation in the mesh refinement algorithm. (Note that, for time-independent test-cases 
such as a Schwarzschild black-hole, the time-interpolation error will effectively not be seen 
due to the small coefficient in front of the corresponding 0(dt 3 ) error term.) 

New insight into the peculiar features of apparent horizons have been revealed in recent 
numerical simulations ll49l and some mathematical machinery has been developed to deal 
with their properties in a rigorous way ll50l IBTl . In particular, a recent theorem due to 
L. Andersson and J. Metzger (L. Andersson, private communication) requires that an outer 
common horizon must already exist if two MOTS come into contact. It is reassuring and 
stimulating that our results are consistent with this theorem. 

5. Conclusions 

Over the last few years substantial progress has been made in developing mathematical 
theorems which establish the well-posedness of the harmonic initial-boundary value problem 
and the stability of its finite difference approximations. We have incorporated some of this 
theory in developing an explicit in time, finite-difference harmonic code and for which we 
have presented a series of tests assessing its validity and accuracy. Such tests range from 
the evolution of isolated black holes to the head-on collision of two black holes and over 
to a binary black hole inspiral and merger. All of them indicate that stable, convergent and 
accurate solutions of the Einstein equations in fully nonlinear regimes can be carried out. 
Furthermore, the merger simulations have also revealed that individual MOTS can touch and 
even intersect. This novel feature in the dynamics of the MOTS was not found before but is 
consistent with theorems on the properties of MOTS. This finding raises new questions about 
the mathematical features of MOTS, which we look forward to exploring. 

Overall, these initial results are an important indication that the AEI harmonic code is 
capable of contributing to the physical understanding of binary black holes currently being 
achieved by numerical simulation. A key assurance of the validity of numerical results for 
binary black holes is their confirmation by a variety of codes based upon different theoretical 
formulations and numerical methods. 

The work reported here also suggests that there is room for a great deal of improvement 
in the numerical methods we have adopted. We have already emphasized the need for 
implementing a constraint-preserving treatment of the outer boundary. In addition, we 
have carried out very little exploration to optimize the use of constraint adjustments, the 
blending of the subluminal and superluminal algorithm, the extrapolation scheme at the 
extraction boundary, mesh refinement, the addition of numerical dissipation, the choice of 
gauge conditions, etc. It is particularly reassuring that the code remains stable even though 
a large amount of error is being generated at the jagged excision boundary. It appears that 
the outflow nature of the boundary, i.e., its spacelike geometry combined with the lightlike 
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direction of all characteristics of the system, leads to an equilibrium between noise generation 
and its flux out of the grid. 

A number of design choices of our code are based on both the experience gained while 
building the Abigel code and on the published work of Pretorius. It is an optimistic sign that 
our present results could be obtained by using what was in a number of cases our initial choice 
for the various code details. We regard this as evidence of the fundamental robustness of the 
harmonic formulation. 
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Appendix A. Blended subluminal-superluminal evolution 

The evolution system (0 consists of coupled quasilinear wave equations whose numerical 
stability is determined by the principle part. By the principle of frozen coefficients l52l . 
the stability analysis can be reduced to a consideration of the wave equation with shift. 
Although finite difference approximations to the wave equation is a well studied problem, the 
complications introduced by a non-zero shift are peculiar to the black hole excision problem. 
This was first recognized in ll53l . where it was suggested that the superluminal shift introduced 
by tracking the excision boundary could be treated by implicit methods. 

Subsequent studies established the stability of explicit finite-difference algorithms, with 
second order accuracy, for the case of superluminal evolution. This was first achieved for the 
ID wave equation with shift 

g u d\<§> + 2d t d x g xt <5> + g xx d 2 x <$> = 0, (A.l) 

in work by Calabrese |39l and Szilagyi et al. BOl . The standard choice of energy for this 
system, 

E W = \Jdx [(-g"(d t $) 2 + fl "i<),A>r , (A.2) 

gives rise to a norm when the evolution direction dt is timelike. In that subluminal case, 
summation by parts (SBP) can be used to establish stability of the semi-discrete approximation 

g u d 2 <S> + 2g xt D d t <f> + g xx D + D-§ = 0, (A.3) 

where D+, D- and Do are, respectively, the standard forward, backward and centered finite 
difference approximations for d x . This ensures that the numerical error is controlled by an 
estimate for the semi-discrete version of the energy norm For most method of lines 

time integrators, e.g. Runge-Kutta, this estimate extends to the fully discretized system. The 
algorithm ( IA.3b has been extended to the 3D subluminal case to give a stable SBP boundary 
treatment E3. 
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However, the algorithm dA.3b is unstable (and cannot be stabilized by Kreiss-Oliger type 
dissipation) when the evolution is superluminal, i.e., when the shift is large enough so that dt 
is spacelike and 

/ xt\2 

g xx = h xx + ^—f- < 0, (A.4) 

gtt 

where h xx > is the inverse to the spatial metric of the t — const Cauchy hypersurfaces. In 
that case, when the energy is no longer a norm, stability can be based upon the positive 
energy associated with the time-like normal n M to the Cauchy hypersurfaces, 

1 :(g tt d t $ + g tx d x $) 2 + h xx dl$ . (A.5) 



= _ / 



gtt' 

As discussed in the ID case [39. 403, the discretization 

xt 

9 U (d t + ^-H ) 2 $ + h xx D+D_<S> = 0. (A.6) 
g tt 

yields a stable second-order accurate superluminal algorithm. Stable superluminal evolution 
algorithms for the 3D case have been given by Motamed et al. l42l . where the global stability 
of a model black hole excision problem is treated. 

Although a stable boundary treatment for the superluminal algorithm dA.6b has been 
proposed BTl . its extended stencil (due to the Dq operator) makes this complicated and an 
SBP boundary version has not yet been formulated. For this reason we use the 3D version of 
the subluminal algorithm ( IA.3I ) in the neighborhood of the outer boundary and blend it to the 
superluminal algorithm ( IA.6b by introducing the vector 

h» = (g tt ,wg it ) (A.7) 

and the evolution variable 

Q = n^$, (A.8) 

where w(x l ) is a spherically symmetric smooth blending function, with w = near the outer 
boundary and w = 1 (so that fit 1 = n^) in the interior. It suffices to discuss the frozen 
coefficient case in which the ID wave equation (IA.U gives rise to the evolution system for Q 
and $, 

g tt d t Q = - (2g xt ~ h x )d x (Q ~ n x d x <S>) - g u g xx d 2 x <S> 

g u d t <S> = Q- n x d x $. (A.9) 

Note that introduction of the auxiliary variable Q, which reduces the system to first-order in 
time, introduces no associated constraints. 

For a second-order accurate approximation, we discretize iA.9i according to 

g u d t Q = -(2g xt - h x )D Q Q + (2g xt - h x )h x D+DS> - g u g xx D + D^ 
g u dt<i> = Q-h x D a <f>. (A. 10) 



In the neighborhood of the outer boundary, this reduces to the subluminal algorithm (IA.3b and 
in the interior where w = 1 it reduces to the superluminal algorithm (IA.6b . The harmonic code 
uses a fourth-order accurate version of (lA.lOb in the interior region. An alternative scheme 
for switching between stable subluminal and superluminal algorithms across the "artificial 
horizon" where dct^g 1 - 7 ) = is given in 
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